@hod_common
pro hod_run,paper=paper,sample=sample,run=run,nburn=nburn,nsamp=nsamp
  COMMON hod_data
  COMMON hod_model
  
if getenv("OMP_NUM_THREADS") eq '' then message,"Set OMP_NUM_THREADS"
splog,"Using ",getenv("OMP_NUM_THREADS")," threads"

; read the observed data
hod_readobs,paper=paper,sample=sample
splog,'read observed data'

parfile=paper+"_"+sample+"_run"+string(run,format='(I02)')+".par"
hod_makeparfile,parfile,rbin_ar=rp_obs,z=meanz_obs ; evaluating model at the observed values of rp
; set the model and parameter file
hod_setmodel,priortype=[1,1,1,1,1,0],$
             priorspec=[ [11.0,15.0],[0.0,1.0],[8.0,15.0],[12.0,15.0],[0.0,1.5],[15.0,0.0] ],$
             parfile=parfile


splog,'set model'


; run MCMC
nstep=10 
if n_elements(nburn) eq 0 then nburn=2000
if n_elements(nsamp) eq 0 then nsamp=10000

if run gt 90 then begin
    ; testing with very small run
    nstep=10
    if n_elements(nburn) eq 0 then     nburn=10
    if n_elements(nsamp) eq 0 then     nsamp=20
endif

hod_mcmc,mean=parmean,error=parerror,samps=samps,lik=lik,/time,/verbose,nstep=nstep,nburn=nburn,nsamp=nsamp,$
  savfile=getenv("HOD_DIR")+'/results/'+paper+'_'+sample+'_run'+string(run,format='(I02)')+'_mcmc.sav'
splog,'finished MCMC'

; save results
savfile=getenv("HOD_DIR")+'/results/'+paper+'_'+sample+'_run'+string(run,format='(I02)')+'.sav'
;save,parfile,nstep,nburn,nsamp,parmean,parerror,samps,lik,filename=savfile
save,filename=savfile,/variables,/comm
splog,'saved results to ',savfile

end
